{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbsphinx": "hidden"
   },
   "outputs": [],
   "source": [
    "import open3d as o3d\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import copy\n",
    "import os\n",
    "import sys\n",
    "\n",
    "# only needed for tutorial, monkey patches visualization\n",
    "sys.path.append('..')\n",
    "import open3d_tutorial as o3dtut\n",
    "# change to True if you want to interact with the visualization windows\n",
    "o3dtut.interactive = not \"CI\" in os.environ"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Point cloud\n",
    "This tutorial demonstrates basic usage of a point cloud.\n",
    "\n",
    "## Visualize point cloud\n",
    "The first part of the tutorial reads a point cloud and visualizes it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Load a ply point cloud, print it, and render it\")\n",
    "ply_point_cloud = o3d.data.PLYPointCloud()\n",
    "pcd = o3d.io.read_point_cloud(ply_point_cloud.path)\n",
    "print(pcd)\n",
    "print(np.asarray(pcd.points))\n",
    "o3d.visualization.draw_geometries([pcd],\n",
    "                                  zoom=0.3412,\n",
    "                                  front=[0.4257, -0.2125, -0.8795],\n",
    "                                  lookat=[2.6172, 2.0475, 1.532],\n",
    "                                  up=[-0.0694, -0.9768, 0.2024])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`read_point_cloud` reads a point cloud from a file. It tries to decode the file based on the extension name. For a list of supported file types, refer to [File IO](file_io.ipynb).\n",
    "\n",
    "`draw_geometries` visualizes the point cloud. Use a mouse/trackpad to see the geometry from different view points.\n",
    "\n",
    "It looks like a dense surface, but it is actually a point cloud rendered as surfels. The GUI supports various keyboard functions. For instance, the `-` key reduces the size of the points (surfels).\n",
    "\n",
    "<div class=\"alert alert-info\">\n",
    "    \n",
    "**Note:** \n",
    "\n",
    "Press the `H` key to print out a complete list of keyboard instructions for the GUI. For more information of the visualization GUI, refer to [Visualization](visualization.ipynb) and [Customized visualization](../visualization/customized_visualization.rst).\n",
    "\n",
    "</div>\n",
    "\n",
    "<div class=\"alert alert-info\">\n",
    "    \n",
    "**Note:** \n",
    "\n",
    "On macOS, the GUI window may not receive keyboard events. In this case, try to launch Python with `pythonw` instead of `python`.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Voxel downsampling\n",
    "Voxel downsampling uses a regular voxel grid to create a uniformly downsampled point cloud from an input point cloud. It is often used as a pre-processing step for many point cloud processing tasks. The algorithm operates in two steps:\n",
    "\n",
    "1. Points are bucketed into voxels.\n",
    "2. Each occupied voxel generates exactly one point by averaging all points inside."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Downsample the point cloud with a voxel of 0.05\")\n",
    "downpcd = pcd.voxel_down_sample(voxel_size=0.05)\n",
    "o3d.visualization.draw_geometries([downpcd],\n",
    "                                  zoom=0.3412,\n",
    "                                  front=[0.4257, -0.2125, -0.8795],\n",
    "                                  lookat=[2.6172, 2.0475, 1.532],\n",
    "                                  up=[-0.0694, -0.9768, 0.2024])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Vertex normal estimation\n",
    "Another basic operation for point cloud is point normal estimation.\n",
    "Press `N` to see point normals. The keys `-` and `+` can be used to control the length of the normal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Recompute the normal of the downsampled point cloud\")\n",
    "downpcd.estimate_normals(\n",
    "    search_param=o3d.geometry.KDTreeSearchParamHybrid(radius=0.1, max_nn=30))\n",
    "o3d.visualization.draw_geometries([downpcd],\n",
    "                                  zoom=0.3412,\n",
    "                                  front=[0.4257, -0.2125, -0.8795],\n",
    "                                  lookat=[2.6172, 2.0475, 1.532],\n",
    "                                  up=[-0.0694, -0.9768, 0.2024],\n",
    "                                  point_show_normal=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`estimate_normals` computes the normal for every point. The function finds adjacent points and calculates the principal axis of the adjacent points using covariance analysis.\n",
    "\n",
    "The function takes an instance of `KDTreeSearchParamHybrid` class as an argument. The two key arguments `radius = 0.1` and `max_nn = 30` specifies search radius and maximum nearest neighbor. It has 10cm of search radius, and only considers up to 30 neighbors to save computation time.\n",
    "\n",
    "<div class=\"alert alert-info\">\n",
    "    \n",
    "**Note:** \n",
    "\n",
    "The covariance analysis algorithm produces two opposite directions as normal candidates. Without knowing the global structure of the geometry, both can be correct. This is known as the normal orientation problem. Open3D tries to orient the normal to align with the original normal if it exists. Otherwise, Open3D does a random guess. Further orientation functions such as `orient_normals_to_align_with_direction` and `orient_normals_towards_camera_location` need to be called if the orientation is a concern.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Access estimated vertex normal\n",
    "Estimated normal vectors can be retrieved from the `normals` variable of `downpcd`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Print a normal vector of the 0th point\")\n",
    "print(downpcd.normals[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To check out other variables, please use `help(downpcd)`. Normal vectors can be transformed as a numpy array using `np.asarray`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Print the normal vectors of the first 10 points\")\n",
    "print(np.asarray(downpcd.normals)[:10, :])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check [Working with NumPy](working_with_numpy.ipynb) for more examples regarding numpy arrays."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Crop point cloud"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Load a polygon volume and use it to crop the original point cloud\")\n",
    "demo_crop_data = o3d.data.DemoCropPointCloud()\n",
    "pcd = o3d.io.read_point_cloud(demo_crop_data.point_cloud_path)\n",
    "vol = o3d.visualization.read_selection_polygon_volume(demo_crop_data.cropped_json_path)\n",
    "chair = vol.crop_point_cloud(pcd)\n",
    "o3d.visualization.draw_geometries([chair],\n",
    "                                  zoom=0.7,\n",
    "                                  front=[0.5439, -0.2333, -0.8060],\n",
    "                                  lookat=[2.4615, 2.1331, 1.338],\n",
    "                                  up=[-0.1781, -0.9708, 0.1608])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`read_selection_polygon_volume` reads a json file that specifies polygon selection area. `vol.crop_point_cloud(pcd)` filters out points. Only the chair remains."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Paint point cloud"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Paint chair\")\n",
    "chair.paint_uniform_color([1, 0.706, 0])\n",
    "o3d.visualization.draw_geometries([chair],\n",
    "                                  zoom=0.7,\n",
    "                                  front=[0.5439, -0.2333, -0.8060],\n",
    "                                  lookat=[2.4615, 2.1331, 1.338],\n",
    "                                  up=[-0.1781, -0.9708, 0.1608])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`paint_uniform_color` paints all the points to a uniform color. The color is in RGB space, [0, 1] range."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Point cloud distance\n",
    "Open3D provides the method `compute_point_cloud_distance` to compute the distance from a source point cloud to a target point cloud. I.e., it computes for each point in the source point cloud the distance to the closest point in the target point cloud.\n",
    "\n",
    "In the example below we use the function to compute the difference between two point clouds. Note that this method could also be used to compute the Chamfer distance between two point clouds."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load data\n",
    "demo_crop_data = o3d.data.DemoCropPointCloud()\n",
    "pcd = o3d.io.read_point_cloud(demo_crop_data.point_cloud_path)\n",
    "vol = o3d.visualization.read_selection_polygon_volume(demo_crop_data.cropped_json_path)\n",
    "chair = vol.crop_point_cloud(pcd)\n",
    "\n",
    "dists = pcd.compute_point_cloud_distance(chair)\n",
    "dists = np.asarray(dists)\n",
    "ind = np.where(dists > 0.01)[0]\n",
    "pcd_without_chair = pcd.select_by_index(ind)\n",
    "o3d.visualization.draw_geometries([pcd_without_chair],\n",
    "                                  zoom=0.3412,\n",
    "                                  front=[0.4257, -0.2125, -0.8795],\n",
    "                                  lookat=[2.6172, 2.0475, 1.532],\n",
    "                                  up=[-0.0694, -0.9768, 0.2024])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bounding volumes\n",
    "The `PointCloud` geometry type has bounding volumes as all other geometry types in Open3D. Currently, Open3D implements an `AxisAlignedBoundingBox` and an `OrientedBoundingBox` that can also be used to crop the geometry."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "aabb = chair.get_axis_aligned_bounding_box()\n",
    "aabb.color = (1, 0, 0)\n",
    "obb = chair.get_oriented_bounding_box()\n",
    "obb.color = (0, 1, 0)\n",
    "o3d.visualization.draw_geometries([chair, aabb, obb],\n",
    "                                  zoom=0.7,\n",
    "                                  front=[0.5439, -0.2333, -0.8060],\n",
    "                                  lookat=[2.4615, 2.1331, 1.338],\n",
    "                                  up=[-0.1781, -0.9708, 0.1608])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Convex hull\n",
    "The convex hull of a point cloud is the smallest convex set that contains all points. Open3D contains the method `compute_convex_hull` that computes the convex hull of a point cloud. The implementation is based on [Qhull](http://www.qhull.org/).\n",
    "\n",
    "In the example code below we first sample a point cloud from a mesh and compute the convex hull that is returned as a triangle mesh. Then, we visualize the convex hull as a red `LineSet`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bunny = o3d.data.BunnyMesh()\n",
    "mesh = o3d.io.read_triangle_mesh(bunny.path)\n",
    "mesh.compute_vertex_normals()\n",
    "\n",
    "pcl = mesh.sample_points_poisson_disk(number_of_points=2000)\n",
    "hull, _ = pcl.compute_convex_hull()\n",
    "hull_ls = o3d.geometry.LineSet.create_from_triangle_mesh(hull)\n",
    "hull_ls.paint_uniform_color((1, 0, 0))\n",
    "o3d.visualization.draw_geometries([pcl, hull_ls])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## DBSCAN clustering\n",
    "Given a point cloud from e.g. a depth sensor we want to group local point cloud clusters together. For this purpose, we can use clustering algorithms. Open3D implements DBSCAN [\\[Ester1996\\]](../reference.html#Ester1996) that is a density based clustering algorithm. The algorithm is implemented in `cluster_dbscan` and requires two parameters: `eps` defines the distance to neighbors in a cluster and `min_points` defines the minimum number of points required to form a cluster. The function returns `labels`, where the label `-1` indicates noise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ply_point_cloud = o3d.data.PLYPointCloud()\n",
    "pcd = o3d.io.read_point_cloud(ply_point_cloud.path)\n",
    "\n",
    "with o3d.utility.VerbosityContextManager(\n",
    "        o3d.utility.VerbosityLevel.Debug) as cm:\n",
    "    labels = np.array(\n",
    "        pcd.cluster_dbscan(eps=0.02, min_points=10, print_progress=True))\n",
    "\n",
    "max_label = labels.max()\n",
    "print(f\"point cloud has {max_label + 1} clusters\")\n",
    "colors = plt.get_cmap(\"tab20\")(labels / (max_label if max_label > 0 else 1))\n",
    "colors[labels < 0] = 0\n",
    "pcd.colors = o3d.utility.Vector3dVector(colors[:, :3])\n",
    "o3d.visualization.draw_geometries([pcd],\n",
    "                                  zoom=0.455,\n",
    "                                  front=[-0.4999, -0.1659, -0.8499],\n",
    "                                  lookat=[2.1813, 2.0619, 2.0999],\n",
    "                                  up=[0.1204, -0.9852, 0.1215])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-info\">\n",
    "    \n",
    "**Note:** \n",
    "\n",
    "This algorithm precomputes all neighbors in the epsilon radius for all points. This can require a lot of memory if the chosen epsilon is too large.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plane segmentation\n",
    "Open3D also supports segmententation of geometric primitives from point clouds using RANSAC. To find the plane with the largest support in the point cloud, we can use `segment_plane`. The method has three arguments: `distance_threshold` defines the maximum distance a point can have to an estimated plane to be considered an inlier, `ransac_n` defines the number of points that are randomly sampled to estimate a plane, and `num_iterations` defines how often a random plane is sampled and verified. The function then returns the plane as $(a,b,c,d)$ such that for each point $(x,y,z)$ on the plane we have $ax + by + cz + d = 0$. The function further returns a list of indices of the inlier points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pcd_point_cloud = o3d.data.PCDPointCloud()\n",
    "pcd = o3d.io.read_point_cloud(pcd_point_cloud.path)\n",
    "\n",
    "plane_model, inliers = pcd.segment_plane(distance_threshold=0.01,\n",
    "                                         ransac_n=3,\n",
    "                                         num_iterations=1000)\n",
    "[a, b, c, d] = plane_model\n",
    "print(f\"Plane equation: {a:.2f}x + {b:.2f}y + {c:.2f}z + {d:.2f} = 0\")\n",
    "\n",
    "inlier_cloud = pcd.select_by_index(inliers)\n",
    "inlier_cloud.paint_uniform_color([1.0, 0, 0])\n",
    "outlier_cloud = pcd.select_by_index(inliers, invert=True)\n",
    "o3d.visualization.draw_geometries([inlier_cloud, outlier_cloud],\n",
    "                                  zoom=0.8,\n",
    "                                  front=[-0.4999, -0.1659, -0.8499],\n",
    "                                  lookat=[2.1813, 2.0619, 2.0999],\n",
    "                                  up=[0.1204, -0.9852, 0.1215])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Planar patch detection\n",
    "In addition to finding the single plane with the largest support, Open3D includes an algorithm which uses a robust statistics-based approach for planar patch detection [\\[ArujoAndOliveira2020\\]](../reference.html#ArujoAndOliveira2020). This algorithm first subdivides the point cloud into smaller chunks (using an octree), then attempts to fit a plane to each chunk. If the plane passes a robust planarity test, then it is accepted. A plane is fitted to a subset of points by taking the median point position and the median point normal and estimating a plane $ax + by + cz + d = 0$. The robust planarity check consists of two main components. First, the distribution of angles between each point normal and the fitted plane normal is found. If the spread of this distribution is too high (i.e., there is too much variance amongst all of the associated point normals), then the plane is rejected. Second, the distribution of the distances from the fitted plane to each point is computed. If the spread of this distribution is too high, as measured using a coplanarity metric (see Fig. 4 of [\\[ArujoAndOliveira2020\\]](../reference.html#ArujoAndOliveira2020)), then the plane is rejected. After an initial set of planes are found, an iterative procedure is used to grow and merge planes into a smaller, stable set of planes. These planes can then be bounded using the 2D convex hull of their associated point sets and planar patches are extracted.\n",
    "\n",
    "To find a list of planar patches in a point cloud, we can use `detect_planar_patches`. This method can take six arguments: `normal_variance_threshold_deg` controls the amount of variance allowed amongst the point normals and takes $60^\\circ$ as its default. Smaller values tends to result in fewer, higher quality planes. `coplanarity_deg` controls the allowed distribution of point distances from the plane and takes $75^\\circ$ as its default. Larger values encourage a tighter distribution of points around the fitted plane. `outlier_ratio` sets the maximum allowable outlier ratio in a fitted planes associated set of points before being rejected and has the default value of `0.75`. `min_plane_edge_length` is used to reject false positives---a planar patch's largest edge much be greater than this value to be considered a true planar patch. If left at `0`, the algorithm defaults to using a value of 1\\% of the point clouds largest dimension. `min_num_points` determines how deep the associated octree becomes and how many points must be present when attempting to fit a plane. If left at `0`, the algorithm defaults to 0.1\\% of the number of points in the point cloud. `search_param` is an instance of `geometry::KDTreeSearchParam` and defaults to `geometry::KDTreeSearchParamKNN`. The `k` nearest neighbors to each point are used when growing and merging planes. Larger values of `k` tend to produce higher quality patches at the expense of computation. This function then returns a list of detected planar patches, represented as `geometry::OrientedBoundingBox` objects, with the third column (i.e., $z$) of `R` indicating the planar patch normal vector. The extent in the $z$ direction is non-zero so that the OrientedBoundingBox contains the points that contribute to the plane detection. The planar patches can be visualized using the `geometry::TriangleMesh::CreateFromOrientedBoundingBox` factory function, using the `scale` parameter to \"flatten\" the bounding box along the normal, i.e., `CreateFromOrientedBoundingBox(obox, scale=[1, 1, 0.0001])`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataset = o3d.data.PCDPointCloud()\n",
    "pcd = o3d.io.read_point_cloud(dataset.path)\n",
    "assert (pcd.has_normals())\n",
    "\n",
    "# using all defaults\n",
    "oboxes = pcd.detect_planar_patches(\n",
    "    normal_variance_threshold_deg=60,\n",
    "    coplanarity_deg=75,\n",
    "    outlier_ratio=0.75,\n",
    "    min_plane_edge_length=0,\n",
    "    min_num_points=0,\n",
    "    search_param=o3d.geometry.KDTreeSearchParamKNN(knn=30))\n",
    "\n",
    "print(\"Detected {} patches\".format(len(oboxes)))\n",
    "\n",
    "geometries = []\n",
    "for obox in oboxes:\n",
    "    mesh = o3d.geometry.TriangleMesh.create_from_oriented_bounding_box(obox, scale=[1, 1, 0.0001])\n",
    "    mesh.paint_uniform_color(obox.color)\n",
    "    geometries.append(mesh)\n",
    "    geometries.append(obox)\n",
    "geometries.append(pcd)\n",
    "\n",
    "o3d.visualization.draw_geometries(geometries,\n",
    "                                  zoom=0.62,\n",
    "                                  front=[0.4361, -0.2632, -0.8605],\n",
    "                                  lookat=[2.4947, 1.7728, 1.5541],\n",
    "                                  up=[-0.1726, -0.9630, 0.2071])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Hidden point removal\n",
    "Imagine you want to render a point cloud from a given view point, but points from the background leak into the foreground because they are not occluded by other points. For this purpose we can apply a hidden point removal algorithm. In Open3D the method by [\\[Katz2007\\]](../reference.html#Katz2007) is implemented that approximates the visibility of a point cloud from a given view without surface reconstruction or normal estimation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Convert mesh to a point cloud and estimate dimensions\")\n",
    "armadillo = o3d.data.ArmadilloMesh()\n",
    "mesh = o3d.io.read_triangle_mesh(armadillo.path)\n",
    "mesh.compute_vertex_normals()\n",
    "\n",
    "pcd = mesh.sample_points_poisson_disk(5000)\n",
    "diameter = np.linalg.norm(\n",
    "    np.asarray(pcd.get_max_bound()) - np.asarray(pcd.get_min_bound()))\n",
    "o3d.visualization.draw_geometries([pcd])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Define parameters used for hidden_point_removal\")\n",
    "camera = [0, 0, diameter]\n",
    "radius = diameter * 100\n",
    "\n",
    "print(\"Get all points that are visible from given view point\")\n",
    "_, pt_map = pcd.hidden_point_removal(camera, radius)\n",
    "\n",
    "print(\"Visualize result\")\n",
    "pcd = pcd.select_by_index(pt_map)\n",
    "o3d.visualization.draw_geometries([pcd])"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Edit Metadata",
  "kernelspec": {
   "display_name": "Python 3.8.13 ('open3d')",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.13"
  },
  "vscode": {
   "interpreter": {
    "hash": "d3ff712f491db3b4ccd20005f84e312b261650c105d83adbc3f83f5ecf60b94b"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
